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Abstract: 

This  report  is  divided  into  three  parts: 

In  the  first  part  we  present  our  ongoing  work  on  the  optimization  of  an  MHD  energy  by-pass  concept. 
Here  we  consider  the  optimization  of  the  power  generator  and  accelerator  components  individually,  and 
are  in  the  process  of  a  simultaneous  optimization  of  an  integrated  generator-combustor-accelerator 
concept  in  a  2-D  sense.  We  have  concentrated  our  efforts  on  developing  an  optimization  scheme  that 
couples  a  flow  solver  (perfect  gas  Euler  and  equilibrium  gas  N-S)  with  a  Poisson  solver  for  the  electric 
field  including  Hall  effects.  The  architecture/algorithm  of  the  optimization  scheme  is  such  that  geometric 
and/or  physical  parameters  can  be  optimized  for  a  given  set  of  free-stream  conditions  and  objective 
function.  The  objective  function  was  MHD  power  extracted  in  the  case  of  a  MHD  generator,  and  thrust  in 
the  case  of  an  accelerator. 

The  second  part  of  this  report  presents  some  ideas  on  how  to  extend  this  development  and  the  associated 
real-gas  MHD  technology  at  HyPerComp  Inc.  into  a  potential  Phase-II.  We  have  developed  a  higher  (4th 
order  and  beyond,)  order  accurate  solver  for  MHD  developed  under  an  AFRL  contract.  We  consider 
possibilities  involving  the  usage  of  this  solver  in  accurate  boundary  layer  calculations  and  plasma  effects 
in  shear  layers  as  a  potential  Phase-II  extension.  The  current  study  of  energy  bypass  concepts  may  itself 
be  extended  into  an  extensive  exploration  of  finite  rate  processes  in  such  systems  coupled  with  an 
efficient  optimization  routine  based  on  adjoint  methods. 

A  masters  thesis  supported  in  part  by  this  contract  on  the  optimization  problem  setup  for  hypersonic  inlets 
to  improve  mass  capture  has  been  completed.  Relevant  portions  of  this  thesis  have  been  appended  to  this 
report  in  the  third  part. 
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1.1  MHD  with  energy  bypass  -  the  associated  optimization  problem 

At  HyPerComp,  we  have  in  the  past  studied  global  flow  control  for  accelerators  and  hypersonic  inlets 
using  MHD.  Further,  there  are  studies  from  the  literature  that  provide  analytical  results  in  MHD  power 
generator  optimization.  Test  data  to  validate  numerical  models  is  available  in  plenty  due  to  the  interest  in 
this  subject  in  the  1960s  and  to  some  small  extent,  even  in  the  present  time.  Due  to  the  high  confidence 
level  in  the  prediction  and  optimization  of  MHD  power  generator/accelerator  design,  as  well  as  some 
recent  AFRL  interest  in  energy  bypass  concepts  (Refs  [3,6]),  we  resolved  to  choose  the  following 
problem  as  a  demonstration  of  Phase-I  progress: 


Hall-type  generator  Hall-type  accelerator 


Figure- 1 :  Energy  bypass  concept  for  hypersonics  with  an  integrated  MHD  generator/accelerator 


An  integrated  MHD  power  generator/accelerator  system  coupled  through  a  combustion  chamber  where  a 
fixed  amount  of  heat  is  added  to  the  system  will  be  optimized  to  obtain  the  best  geometric  parameters. 
Applied  magnetic  field  will  be  fixed,  and  a  Hall-type  system  will  be  used  to  obtain  best  performance  from 
the  device.  An  equilibrium  air  Navier-Stokes  solver  will  be  coupled  with  a  Poisson  solver  which  will 
compute  the  electric  field.  Due  to  the  need  for  fine  segmentation  in  the  electrodes  to  overcome  the  loss  of 
efficiency  due  to  the  Hall  effect,  a  special  set  of  boundary  conditions  will  be  used  at  electrode  walls. 
Electrical  conductivity  produced  by  seed  substances  as  well  as  non-equilibrium  e-beam  ionization  can  be 
studied  using  our  codes.  However,  we  will  restrict  Phase-I  work  to  include  only  equilibrium  reacting  air 
with  alkali  seed.  The  team  from  UT-Arlington  will  conclude  their  investigation  on  adjoint  techniques 
initiated  in  the  earlier  progress  report.  Further,  they  will  perform  trial  studies  on  hypersonic  inlets  with 
MHD  wherein  they  will  attempt  to  optimize  the  mass  capture  at  off-design  conditions  using  an  applied 
magnetic  field  in  the  plane  of  the  flow. 


The  success  of  Phase-I  research  will  demonstrate  the  use  of  optimization  techniques  in  the  design  and 
estimation  of  the  potential  performance  of  energy  bypass  devices  being  proposed  for  hypersonic  flight.  It 
will  also  produce  a  problem  formulation  which  may  be  extended  in  Phase-11  to  study  alternate  “local” 
MHD  flow  control  applications  that  are  beyond  the  scope  of  phase-I  study. 
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1.2  MHD  power  generator  /  accelerator 


Segmented  electrodes 


Figure  7:  Schematic  for  MHD  power  generator  optimization  study 


The  theory  of  MHD  power  generators  and  that  of  flow  accelerators  is  very  similar  in  nature  and  is  in  a 
fairly  mature  stage,  having  been  an  active  area  of  research  since  the  1960s  (the  ideas  date  back  to 
Faraday!).  Interest  in  optimizing  their  design  came  from  worldwide  efforts  to  build  them  and  the 
conviction  that  their  efficiency  would  be  high  compared  to  the-then  available  power  producing  methods 
due  to  the  high  temperatures  of  the  working  medium  and  the  ability  to  dispense  with  expensive  rotating 
machinery  such  as  in  turbines.  Starting  from  the  work  of  Neuringer  [11],  several  analytical  studies  (1,15) 
have  been  performed  to  study  the  optimization  problem.  A  recent  study  attempted  to  study  the  same 
problem  (ref.  [5]),  quite  characteristically  unaware  of  the  wealth  of  prior  technical  literature  in  the  subject. 

The  recent  development  of  advanced  modeling  ability  in  the  area  of  MHD  and  plasma  physics  (e.g., 
Gaitonde  [3]  among  others,)  has  enabled  researchers  to  use  CFD  to  assess  the  possibility  of  in-flight 
MHD  power  generation  and  flow  acceleration,  particularly  for  energy  bypass  concepts  as  discussed  in 
refs  [3,6].  We  propose  now,  to  use  the  studies  of  refs  [3,6]  as  a  baseline  and  perform  a  formal  design 
optimization  study,  albeit  2-D  in  phase-I,  to  assess  the  sensitivity  of  its  performance  to  various  design 
parameters. 


In  this  report,  we  merely  present  an  overview  of  the  physics  and  some  sample  results  that  validate  the 
modeling  of  Hall-type  power  generators  using  MHD  codes  at  HyPerComp. 


1.2.1  Ohm’s  law  relevant  to  Hall  generators  and  accelerators 

The  appropriate  form  of  Ohm’s  law  relevant  to  MHD  power  generation  or  flow  acceleration  must  include 
the  Hall  effect  and  may  be  then  written  as: 


J  =  o(E  +  V  xB) - (JxB) 

B 

The  electric  field  (applied  or  induced,)  must  be  divergence  and  curl-free  and  can  be  derived  from  a  scalar 
potential.  It  is  also  at  times  customary  to  define  a  total  electric  field  as  the  field  observed  by  a  particle 
traversing  with  the  local  velocity  of  the  medium.  These  statements  can  be  expressed  mathematically  as: 

E  =  -V<p,  E*  -  E  +  V  x  B 

We  assume  a  two  dimensional  flow,  in  which  velocity  and  electric  field  components  only  have  x  and  y 
components.  At  each  domain  boundary  or  computational  cell  face,  a  normal  vector  and  an  accompanying 
tangent  vector  may  be  defined.  We  use  the  following  notation  for  these  vector  quantities: 

E'  =(E*x,E'y),  V  =  (u,v),  h  =  (nx,ny\  x  =  (rx,ty)  =  (-  ny,nx) 
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The  final  form  of  Ohm’s  law  that  is  used  in  our  computations  is  written  as  follows: 


u 

1 

ha 

h  * 

ji 

_ l 

-<j 

(<PX  ~  Bv)-  fi(<py  +  Buf 

1  +  J32 

E)  +  PE'X 

~  \  +  p2 

[<Py  +  Bu)+  p(tpx  -  Bv) 

The  divergence  of  the  current  density  J  derived  above  must  identically  vanish  for  an  electrically  neutral 
fluid.  This  results  in  a  Poisson  type  equation  for  the  electric  potential  which  must  be  solved  for  an 
appropriate  set  of  boundary  conditions.  We  present  some  examples  of  such  boundary  conditions. 

Electrode  surface:  An  electrode  is  a  constant  potential  surface.  To  this  extent,  one  may  clamp  the  potential 
of  any  one  electrode  in  the  domain  to  be  the  ground  potential  zero.  However,  any  other  electrodes  present, 
must  inherit  a  potential  based  on  the  exact  layout  of  an  intemal/extemal  circuit  connecting  these 
electrodes.  In  the  case  of  an  accelerator,  a  fixed  potential  difference  may  be  applied  across  pairs  of 
electrodes.  While  the  potential  of  one  electrode  in  a  pair  may  be  computed  from  an  appropriate  BC,  the 
other  is  obtained  by  adding  the  voltage  of  the  battery.  In  a  generator,  a  load  may  be  connected  across  a 
pair  of  electrodes  and  the  potential  difference  may  be  inferred  from  the  current  circuited  externally 
through  the  electrode  pair. 

Insulator:  An  insulating  wall  has  zero  current  flowing  normally  into  it.  This  can  be  used  as  a  BC. 

Infinitely  segmented  electrodes:  The  Hall  effects  diminishes  the  performance  of  MHD  devices  by 
producing  current  components  normal  to  the  principal  direction  of  current  production  and  diminishing  the 
magnitude  of  the  total  current  by  a  factor.  Segmentation  of  the  electrodes  can  somewhat  alleviate  this 
effect  and  improve  performance.  In  a  practical  CFD  code,  modeling  a  very  finely  segmented  set  of 
electrodes  is  frequently  infeasible  due  to  the  large  demands  on  the  mesh  resolution.  This  issue  can  be 
resolved  by  the  assumption  of  infinite  segmentation,  as  discussed  in  Oliver  and  Mitchner  [12],  The 
boundary  condition  to  be  used  in  such  a  case  is  that  the  current  tangential  to  the  region  with  infinitely 
segmented  electrodes  is  zero.  Further,  these  electrodes  may  be  connected  externally  in  order  to  reduce  the 
transverse  potential  gradient  in  a  “Linear  Hall  Generator”.  In  such  a  situation,  the  potential  of  one  row  of 
electrodes  is  let  to  float  subject  to  vanishing  tangential  current.  The  corresponding  set  on  the  opposite  face 
are  given  a  potential  which  is  equal  to  it  or  determined  from  external  loads  connected  between  them. 

1.2.2  Mathematical  expressions  for  current  BCs: 

As  mentioned  earlier,  we  use  two  types  of  BCs  on  current  related  to  its  component  normal  or  tangential  to 
a  wall.  Using  the  notation  from  before,  we  may  write  these  components  as: 


J  •  ft  =  —^-((pxnx  -  Bvnx  -  Pq>  nx  -  (3Bunx  +  (pyny  +  Buny  +  ptpji  -  0Bvn  ) 


1  +  J3 


=  ($»*("*  +  fa,)+  <Py(ny  -  fax)-  B(v  ■  t)-(3B{V  • «)) 


J  T  =  - ~^p{(p\fax  -  ny)+  (py{nx  +  pny)+  B(V  ■  n)- /3B{V  ■  f)) 

In  the  case  when  these  components  are  to  vanish,  we  can  obtain  convenient  expressions  for  the  gradient  of 
the  electric  potential  by  rearranging  terms,  as  below: 
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J-T  =  0=><p 
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The  case  of  a  linear  Hall  generator  has  been  studied  often  in  the  past.  Here,  an  infinite  sequence  of 
electrode  pairs  is  short-circuited  from  outside.  This  makes  the  potential  difference  in  the  y-direction  zero: 

Ey=  0 

Current  components  can  then  be  obtained  by  the  somewhat  simplified  expressions: 

Jx  =  — (Ex  +  UJ3B),  J  =  — ^  (J3EX  -  uB) 

Further  simplification  may  result  when  studying  the  open-circuit  condition  wherein  the  current  in  the  x- 
direction  is  zero  in  a  one-dimensional  sense.  All  field  components  can  be  derived  analytically  then. 

1.2.3  Faraday-type  MHD  channel  with  Hall  effect 

We  first  present  an  MHD  channel  with  a  Hall  parameter  /?=  1  and  electrical  conductivity  =  1  (non- 
dimensional  quantities).  A  Faraday  type  configuration  is  considered,  where  flat  plate  electrodes  are  used 
to  cover  a  significant  portion  of  the  channel.  These  electrodes  are  held  at  a  constant  potential  difference. 
Electric  potential  contours  and  current  lines  for  such  a  flow  are  shown  in  Fig.  [8]  below. 


1  10  It 

Figure  8:  Electric  potential  contours  (left)  and  current  lines  (right) 


A  characteristic  feature  of  such  a  flow  is  the  tendency  of  current  lines  to  bunch  upstream  of  the  cathode 
and  downstream  of  the  anode  in  the  presence  of  the  Hall  effect.  The  strong  concentration  of  field 
gradients  in  these  regions  can  cause  damage  to  the  electrodes,  and  is  one  of  the  perils  of  segmenting 
electrodes  in  MHD  channels. 

While  the  above  describes  the  current  features  at  the  electrode-insulator  boundary,  the  applied  magnetic 
field  has  a  very  direct  bearing  on  the  induced  current,  as  is  shown  in  the  following  cases.  The  applied 
magnetic  field  is  ramped  up  and  held  constant  through  most  of  the  channel.  All  channel  walls  are  assumed 
to  be  insulating  here.  The  magnetic  field  is  then  evenly  ramped  down  to  zero.  The  regions  of  B-field 
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gradients  show  the  presence  of  eddy  currents,  as  seen  in  Fig.  [9],  As  the  Hall  parameter  increases,  these 
features  develop  asymmetries,  and  can  eventually  form  very  strong  current  components  near  the  wall. 


Figure  9:  Current  lines  (above)  and  potential  contours  (below)  for  a  channel  flow  with  B-field  gradients. 
Hall  parameter  is  0,  1  and  2  in  the  3  images,  in  sequence. 


1.2.4  Infinitely  segmented  linear  Hall  generator 


A  very  critical  part  of  our  approach  to  MHD/plasma  flows  is  to  model  only  verifiable  physical 
phenomena  for  which  there  are  adequate  numerical  models  in  the  literature.  The  study  of  linear  Hall 
generators  is  virtually  a  classical  problem  in  MHD  power  generator  studies  and  serves  as  an  excellent  test 
case  in  checking  an  MHD  Poisson  solver  routine  for  boundary  conditions  and  internal  consistencies.  We 
use  the  relations  presented  earlier  for  quasi- ID  analysis  to  cross  check  the  results  from  a  2-D  calculation. 
The  flow  parameters  are  selected  as  follows: 

w  =  3.55,  v  =  0,  <7  =  1,  /?  =  1,  Bmm  =1 
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Figure  10:  B-field  and  electric  potential  distribution  for  the  segmented  electrode  generator 
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We  present  in  Fig  [11]  the  results  from  this  calculation  in  the  region  of  the  Hall  generator.  Using  the 
parameters  listed  earlier,  and  the  relations  in  section  [4.2]  for  a  linear  Hall  generator,  we  must  obtain  a  y- 
direction  current  density  of  0.5(-uB  +  /?E<).  The  component  Ex  has  been  computed  to  be  about  -2.931. 
The  value  of  uB  is  about  -3.55.  This  gives  a  current  density  jy  of  about  -3.24  which  matches  closely  with 
the  mean  value  shown  in  Fig.  [11]. 


x 


Further,  it  may  be  observed  that  the  electric  potential  gradient  in  the  y-direction  is  virtually  zero,  while  a 
Hall  voltage  develops  along  the  length  of  the  channel.  All  these  parameters  are  observed  to  be  in  good 
agreement  with  quasi- ID  estimates,  thus  validating  the  MHD  portion  of  this  modeling  effort. 
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1.3  Shape  Optimization 


In  this  study  we  concentrated  on  optimizing  a  single  geometric  parameter,  in  order  to  demonstrate  the 
feasibility  of  our  approach.  This  study  was  also  intended  to  identify  issues  such  as  convergence  criterion 
for  the  objective  function,  effect  of  grid-size  on  the  search  direction  in  the  optimization  routine  and 
computing  time  for  an  optimization  cycle.  It  was  also  aimed  at  identifying  limits  of  physical  parameters 
(electrical  conductivity/interaction  parameter)  on  the  stability  and  convergence  of  the  overall  scheme. 
Optimization  results  shown  here  are  for  inviscid  flows  (Euler  solutions).  The  results  from  these  solutions 
would  be  good  indicators  to  assess  grid  size  and  operating  range  (for  design  variables)  for  viscous  flow 
simulations. 

Figure  1 2  shows  a  generator/accelerator  formed  using  a  double  cubic  spline.  The  design  parameters  and 
their  notations  are  explained  in  Table  I.  Variation  of  these  design  parameters  can  be  used  to  optimize  the 
shape  of  the  accelerator/generator.  In  the  results  presented  here,  a  single  design  parameter,  (y’0  or  y’j)  has 
been  varied  in  the  optimization  process.  The  constraints  on  the  design  variables  were  that  the  slope  be 
positive. 


/ 


Table  I 


Parameter 

Symbo 

1 

Inlet  Diameter 

Yo 

Outlet  Diameter 

y  / 

Location  of  end-points  of 
spline 

xc 

Slope  at  inlet 

V’n 

Slope  at  outlet 

y’/ 

Slope  at  xc 

X’ 

A  o 

Length  of 

generator/accelerator 

i 

Figure  12:  Schematic  of  geometry 

Figure  13  shows  the  shape  optimization  for  maximizing  thrust  when  the  design  variable  is  y ’/ .  All  other 
design  variables  shown  in  Table  I  were  kept  constant  (see  Table  II).  The  initial  value  of  y ’/  is  10  °.  The 
upper  bound  of  y  ’/  was  set  at  40  °.  The  optimization  process  shows  that  increasing  y ’/  increases  thrust. 
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Figure  13:  Shape  optimization  process,  showing  dependence  of  optimized  shape  on  axial  grids. 
Table  II.  Parameters  used  in  the  shape  optimization  study 


Parameter 

Inlet  Diameter  (y0) 

EEXSSlilfii 

Outlet  Diameter  (y /) 

liEKWi 

Location  of  end-points  of  spline  (xc) 

0.0 

0.0 

10.0 

Slope  at  xc 

0.0 

Length  of  generator/accelerator  (/) 

An  important  consideration  in  shape  optimization  is  determining  the  grid  size  so  as  not  to  influence  the 
search  direction  in  the  optimization  process.  Table  III  shows  the  effect  of  axial  grids  on  the  optimization 
process.  The  results  show  converged  values  of  the  objective  function  (thrust)  for  y/  =  10  degrees  (red 
line  in  Figure  2).  It  is  seen  that  when  the  number  of  axial  grid  points  is  below  200,  the  optimization 
search  direction  seeks  a  different  local  maxima.  A  particular  geometry  (corresponding  to  a  given  y/’)  was 
considered  to  be  converged  (for  all  grids)  if  the  variation  in  the  objective  function  was  less  than  0.5% 
after  an  interval  of  500  time-steps.  Thus  a  fully  converged  objective  function  on  a  coarse  grid  might 
indicate  a  different  optimized  solution  as  compared  to  a  fine  grid.  In  this  particular  case,  the  coarser  grid 
indicates  that  the  thrust  is  maximum  if  the  channel  radius  is  monotonically  increasing  and  the  slope  at  the 
outlet  is  near  zero.  The  fine  grid  solutions  (200  X  60,  400  X  60  and  800  X  60)  indicate  however,  that  the 
thrust  can  be  optimized  if  there  is  a  constriction  in  the  channel  between  the  inlet  and  the  outlet,  thus 
implying  that  higher  values  of  yf,  yield  higher  thrust.  These  considerations  are  important  in  deciding  the 
minimum  number  of  grid  points  to  ensure  unambiguous  designs.  Minimizing  the  number  of  grid  points  is 
of  great  importance  in  coupled  MHD-viscous  flow  simulations  as  this  is  a  computationally  intensive 
process. 
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Table  III:  Variation  of  Thrust  at  the  exit  plane  with  axial  grid  size 


Simulation  # 

Thrust  in  KJ 

1 

100X  60 

0.1945442 

2 

200  X  60 

0.2582707 

3 

400  X  60 

0.2924424 

4 

800  X  60 

0.3096796 

From  Table  III  it  is  clear  that  the  thrust  varies  less  than  4%  when  the  number  of  grid  points  in  the  axial 
direction  is  increased  from  400  to  800.  Hence  this  grid  could  be  used  to  conduct  further  optimization 
studies.  It  was  found  that  use  of  non-uniform  grids  in  the  radial  direction  did  not  change  the  value  of 
thrust  obtained  using  uniform  grids. 

1.3.1  Shape  optimization  of  a  MHD-accelerator  (Euler  solutions) 

Based  on  the  results  in  Table  III,  a  400X60  grid  was  used  to  study  the  shape  optimization  process  with 
and  without  MHD  for  conditions  shown  in  Table  IV.  These  operating  conditions  (temperature/pressure 
and  Mach  number)  are  representative  for  operation  of  a  MHD-bypass  accelerator  where  seeding  can  be 
used  to  create  the  required  ionization  levels  for  MHD  effects. 


Table  IV:  Operating  conditions  for  MHD  accelerator 


Inlet  Pressure 


Inlet  Density 


Inlet  Temperature 


Inlet  Mach  Number 


Electrical  conductivi 


Max  B-field 


Hall  Parameter 


Applied  Potential 
Difference 


Interaction  parameter 


0.1  atm 


1 

k 


2500  K 


3 


1  mho/m 


1  T 


Figure  3  shows  the  variation  of  thrust  with  y)  with  and  without  MHD  effects.  It  is  seen  that  the  search 
direction  with  and  without  MHD  is  the  same.  Optimized  thrust  without  MHD  is  less  than  that  obtained 
with  MHD.  The  electrical  conductivity  was  kept  constant  at  1 .0  mho/m  in  the  flow  field.  The  magnetic 
field  (B)  is  as  shown  in  Figure  14. 

B 


1.0 

,  / 

7 

7 

\  , 

0.15 

0.3 

0.7 

0.85 

Figure  14:  Variation  of  Magnetic  field  along  axis  of  MHD  accelerator 


Final  report  for  Phase-I  STTR  Contract  #  FA9550-04-C-01 17,  “MHD  Design  optimization”  Page  11  of  30 

ij  yPerComp,  Inc. 


Laadar  In  High  Farfonnaoca  Computing 


The  maximum  magnetic  field  is  maintained  in  the  middle  of  the  accelerator  in  order  to  prevent  MHD 
effects  from  creating  instabilities  at  the  inlet  and  outlet.  The  magnetic  field  is  along  the  positive  Z-axis 
(out  of  the  plane  of  the  paper).  Figures  1 5  (a)  and  (b)  show  the  Mach  number  contours  with  and  without 
MHD  effects.  For  a  channel  with  an  area-ratio  of  2  and  an  inlet  Mach  number  =  3.0,  the  isentropic  value 
of  the  exit  plane  Mach  number  is  3.75.  Figure  16  shows  that  the  axial  variation  of  the  Mach  number 
along  the  centerline  of  the  accelerator.  It  is  seen  that  the  exit  plane  Mach  number  is  about  3.5  without 
MHD  and  3.8  with  MHD.  The  increase  in  exit  plane  Mach  number  is  small  on  account  of  the  small  value 
of  the  interaction  parameter.  Simulations  are  currently  being  run  for  higher  values  of  interaction 
parameter. 
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Figure  1 5  showing  Mach  number  contours  without  and  with  MHD 


Figure  16:  Axial  variation  of  Mach  number  along  the  centerline. 

Figure  17  shows  contours  of  electric  potential  (9)  along  with  the  current  lines  in  the  channel.  The 
presence  of  the  Hall  effect  is  depicted  in  the  field  lines.  It  is  seen  that  the  current  flows  in  the  positive  Y- 
direction  thus  producing  a  JXB  force  in  the  positive  X-direction  (along  the  flow)  leading  to  flow 
acceleration. 
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Figure  17:  Contours  of  Electric  Potential  with  current  lines 
1.3.2  Shape  Optimization  of  an  Infinitely  segmented  Hall  generator  generator  (Euler  Solutions) 

Motivation  for  simulating  a  Hall  generator: 

In  a  typical  hypersonic  vehicle  operating  at  Mach  numbers  below  10,  air  entering  the  MHD  channel  is  at 
static  pressures  between  0.01  to  0.1  atm  and  temperatures  around  a  few  hundred  Kelvin.  The  electrical 
conductivity  of  air  under  these  conditions  is  very  low  (around  0.1  mho/m  or  less).  Hence  external  on¬ 
board  seeding  or  non-equilibrium  methods  are  needed  to  generate  and  sustain  electron  number  densities 
required  to  obtain  electrical  conductivities  on  the  order  of  10-30  mho/m.  For  a  meaningful  operation, 
power  cost  of 

ionization  should  be  much  lower  than  the  extracted  electric  power.  This  requirement  suggests  the  use  of 
strong  magnetic  fields  (7-10  T)  in  order  to  generate  an  interaction  parameter  high  enough  for  appreciable 
MHD  effects.  Low  pressures  and  high  magnetic  fields  lead  to  high  Hall  parameters,  making  it  difficult  to 
operate  a  MHD  channel  as  a  Faraday  generator.  A  Faraday  channel,  allowing  a  longitudinal  Hall  current 
to  flow  would  sharply  reduce  performance.  The  Hall  current  could  be  prevented  by  segmenting 
electrodes,  but  it  this  could  still  lead  to  engineering  difficulties  such  as  arcing  and  also  deteriorate  the 
performance  of  the  Faraday  channel. 

Given  these  problems  with  a  Faraday  generator,  it  would  be  beneficial  to  use  Hall  configuration,  where 
opposing  pairs  of  segmented  electrodes  are  short-circuited,  and  the  longitudinal  Hall  current  is  collected. 
A  Hall  generator  also  seems  like  a  natural  choice,  given  the  large  Hall  parameters  characterizing  a  MHD- 
based  hypersonic  vehicle.  Hall  generators  also  promise  to  attenuate  the  electrode  arcing  problem  on 
account  of  the  large  longitudinal  current. 

Results: 

In  practical  MHD  generators,  keeping  the  surface  area  (volume)  of  a  generator  at  a  minimum  is  important 
from  the  standpoint  of  minimizing  heat  losses.  A  sample  optimization  problem  based  on  this  requirement, 
wherein  the  slope  of  the  channel  at  the  inlet  was  a  design  parameter  was  studied.  Based  on  this  angle,  the 
length  of  the  Hall  generator  (region  where  the  top  and  bottom  electrodes  were  shorted)  was  varied,  so  as 
to  keep  the  total  area  of  the  generator  constant.  A  higher  slope  at  the  inlet  would  just  lead  to  a  shorter 
axial  length  of  the  generator  and  vice-versa. 
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Figure  18  (a)  shows  the  contours  of  electric  field  with  current  lines,  whereas  Figure  18  (b)  shows  the 
variation  in  inlet  angle  with  optimization  cycles  to  extract  maximum  power.  The  results  indicate  that  a 
smaller  inlet  angle  (and  hence  a  longer  MHD  generator  length)  would  yield  maximum  power  while 


keeping  the  area  fixed. 


Figure  8(b).  Iterative  optimization  of  inlet  angle 


1.3.3  Viscous  simulation  of  MHD  Accelerator/generator: 

The  above  results  were  obtained  using  an  Euler  solver  for  simulating  the  flow.  These  simulations  are 
important  in  order  to  obtain  a  good  estimate  of  the  axial  grids  required  and  the  upper  and  lower  bounds 
for  various  design  variables  of  interest  to  the  problem.  Viscous  flow  simulations  of  these  high-speed 
flows  require  considerably  longer  time  to  converge,  as  it  is  necessary  to  resolve  the  thin  boundary  layers 
in  these  flows.  A  considerable  amount  of  computing  time  can  be  saved  if  one  were  to  conduct  these  N-S 
simulations  in  the  neighborhood  of  the  actual  optimum  solution. 


As  a  part  of  our  Phase  I  effort,  we  have  developed  a  optimization  scheme  for  viscous  similar  to  the  one 
with  Euler  solutions.  This  scheme  is  coupled  with  a  CEA  code  which  can  determine  the  equilibrium 
chemical  compositions  for  these  MHD  flow  and  compute  realistic,  spatially  varying  electrical 
conductivities  and  Hall  parameters.  The  Poisson  solver  can  thus  simulate  a  realistic  MHD 
generator/accelerator  and  thus  be  of  immense  use  as  a  practical  design  tool. 


Figure  19  shows  the  electric  potential  for  an  accelerator  operating  under  the  same  flow  conditions  as 
mentioned  above.  The  flow  comprises  of  air  seeded  with  about  0.02%  of  KOH  at  2500K  at  the  inlet,  thus 
maintaining  an  electrical  conductivity  of  about  10  mho/m. 
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Figure  19:  Contours  of  Electric  field  in  N-S  simulations. 


1.3.4  Shape  optimization  of  an  MHD  accelerator 

We  initiated  a  study  of  seeded  air  MHD  accelerators  partly  because  the  conditions  therein  are  close  to 
equilibrium  and  an  equilibrium  property  code  such  as  CEA  from  NASA  Lewis  can  be  effectively  used  to 
model  gas  properties  in  such  cases.  Electrical  conductivity  is  obtained  from  collisional  cross  section  data 
and  the  local  values  of  species  concentration  of  the  ionizing  medium  (in  this  case,  air  seeded  with  KOH). 
Air  seeded  in  this  manner  develops  a  fair  amount  of  conductivity  at  temperatures  in  the  vicinity  of  3000K. 

We  pose  the  following  optimization  problem: 

Maximize  the  thrust  produced  by  a  nozzle: 

T=  \Pe~  Pa  +  P“2  & 

exit 


with  a  nozzle  length  of  lm,  inlet  width  of  0.1  m  and  exit  width  of  0.5  m,  by  optimizing  the  contour  of  the 
nozzle  surface  defined  by  two  con-joined  segments  shown  in  fig.  [20].  The  first  is  a  circular  arc  segment 
of  a  radius  of  curvature  0.15  m,  and  the  second  is  a  cubic  curve  connected  to  the  first  at  an  angular 
location  0  and  having  a  y-coordinate  of  yl  at  the  exit  plane.  The  design  variables,  therefore,  are  0  and 

yi- 


Figure  20  :  Nozzle  geometry  for 
design  optimization 
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A  preliminary  design  study  without  the  effect  of  MHD  was  made  for  such  a  nozzle  geometry.  Fig  [21] 
shows  a  sample  convergence  history  of  the  objective  function  and  the  arc  angle  as  the  optimization 
progresses.  Nozzle  inflow  Mach  number  was  selected  to  be  3  and  an  exit  plane  ambient  pressure  was 
selected  to  be  the  isentropic  value  for  a  “mean”  area  ratio  of  5. 


0.65 


0j6 
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0.45 


Figure  21:  Nozzle  thrust  optimization  with  the  arc-cubic  geometry  at  inflow  Mach  number  =  3 

A  central  segment  of  such  a  nozzle  was  then  powered  by  electrodes  and  a  linear  Hall-type  generator  was 
modeled.  Using  the  “mean”  design  condition  referred  to  above  (based  on  isentropic  flow),  we  attempt  to 
obtain  the  load  parameter  K  =  Ey/uB  of  about  1.2.  We  select  gas  parameters  such  that  the  mean 
conductivity  is  about  2  mho/m.  A  magnetic  field  of  1  Tesla  perpendicular  to  the  plane  of  the  flow  was 
imposed  in  the  powered  region  with  linear  taper  on  either  side.  The  value  of  the  applied  electric  field  was 
3600  V/m,  perhaps  corresponding  to  batteries  of  about  360  -  450  Volts  connected  across  the  channel.  We 
present  results  from  this  design  problem  below. 
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Figure  22:  v-component  of  velocity  with  MHD  (left)  and  without  MHD  (right) 
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Figure  23:  Optimization  history  showing  the  arc  angle  variation  with  thrust 


Figure  24:  Exit  plane  velocity 
profile  with  and  without  MHD. 
Note  the  sharper  drop  in  the 
velocity  near  the  wall  with 
MHD. 
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Figure  25:  Mach  number  contours  -  note  that  for  these  conditions,  exit  Mach  number  in  the  MHD  case  is 
actually  slightly  lower  than  without  MHD,  while  the  velocity  is  actually  higher. 


1.4  An  integrated  MHD  generator-accelerator  model 

We  made  a  preliminary  study  of  the  coupled  MHD  generator-accelerator  configuration.  Two  casese  were 
studied,  the  first  being  a  flat  duct  with  specialized  areas  for  power  generation,  combustion  and  flow 
acceleration.  The  second  case  was  that  of  a  piecewise  linear  geometry. 


Case-I:  Straight  duct. 


Lengths: 

Generator  =  2.0  m 
Combustor  =  2.4  m 
Accelerator  =  3.3  m 
Total  length  =  7.7  m 
Width  =  0.5  m 
MHD  setup: 

Sigma  (average)  =  15.0 
Beta  =  1 .0 

Accelerator  PD  =  2000  V 
B-field  =  1  T 
Flow  conditions: 

Inlet  Mach  No  =  3.0 
Inlet  Pressure  =  0. 1  atm 
Inlet  Temp  =  2500  K 
Inlet  Velocity  ~  3000  m/s 


Figure  26:  Contours  of  electric  potential  in  the  straight 
duct  accelerator 


The  above  conditions  result  in  an  interaction  parameter  of  about  0.4.  Figure  above  shows  contours  of 
electric  potential  in  such  a  case.  An  effective  power  generated  by  the  linear  Hall  generator  was  0. 1 2  MW, 
while  the  power  used  in  the  accelerator  was  77  MW.  Clearly,  these  are  non-optimal  and  there  is  a  very 
important  role  in  parameter  selection  for  even  so  simple  a  problem.  If  the  system  is  to  produce  a  net 
positive  power,  geometric  as  well  as  electrical  parameters  need  to  be  selected  correctly. 
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Case-II:  Piecewise  linearly  varying  duct. 


Lengths: 

Generator  =  2.0  m 
Combustor  =  2.4  m 
Accelerator  =  3.3  m 
Total  length  =  7.7  m 

Widths: 

Inlet/Outlet  =  0.6  m 
Combustor  width  =  0.3 
MHD  parameters: 

Sigma  =15.0 
Beta  =  1 .0 

Accelerator  PD  =  2000  V 
B-field  =  1  T 

Interaction  parameter  ~0.4 
Flow  parameters: 

Inlet  Mach  No  =  3.0 
Inlet  Pressure  =  0. 1  atm 
Inlet  Temp  =  2500  K 
Inlet  Velocity  ~  3000  m/s 


There  is  a  vast  potential  for  further  investigations  of  this  problem.  Among  the  most  important,  from  a 
systems  perspective,  would  be  the  electrical  coupling  between  the  power  generator  and  accelerator  and 
how  such  a  system  may  operate  in  a  self-sufficient  mode  at  a  high  efficiency  for  hypersonic/scramjet  type 
applications.  It  may  be  noted  that  we  have  again  used  a  seeded-equilibrium  air  type  approach  here  for  the 
sake  of  convenience.  The  introduction  of  non-equilibrium  ionization  would  require  the  use  of  more 
complex  numerical  tools  that  are  presently  beyond  the  scope  of  this  investigation,  but  would  make  an 
interesting  and  important  research  topic  for  a  Phase-II  project. 
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1.5  Summary  of  Results: 

As  a  part  of  our  Phase  I  research,  we  have  concentrated  on  developing  the  tools  and  strategy  to  optimize 
MHD-based  accelerator/generators.  We  have  demonstrated  the  ability  to  simulate  an  MHD 
generator/accelerator  and  optimize  their  shape  for  a  given  design  objective  (maximizing  power  or  thrust). 
We  have  developed  the  capability  to  include  important  physics  (viscous  effects,  spatially  vaiying 
properties  such  as  electrical  conductivities  and  Hall  parameter)  to  simulate  and  optimize  realistic 
operating  conditions  in  a  MHD  by-pass  concept. 
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Plans  for  extension 


Here,  we  present  some  possibilities  for  the  extension  of  this  work  into  a  second  phase.  By  way  of 
applications,  we  wish  to  consider  either  the  bypass  engine  concept  studied  earlier  in  the  exploratory  level 
in  Phase-I,  or  the  localized  flow  control  concepts  that  are  attractive. 
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2.1  Development  of  an  optimal  MHD-energy  bypass  propulsion  system  design 

The  optimal  design  of  an  MHD-energy  bypass  design  for  hypersonic  vehicles  presented  here  can  be 
extended  into  Phase-II  research.  The  design  space  will  indeed  be  extremely  complex,  due  to  the  vast 
diversity  in  time  and  length  scales  in  the  physics  encountered.  We  present  here  a  sampling  of  design 
variables  that  may  potentially  be  chosen  to  model  this  problem. 

The  optimization  problem  itself  can  be  one  of  the  following: 

1)  Optimization  of  the  total  extracted  power 

2)  Optimization  of  efficiency  at  a  given  power/length  (minimization  of  losses) 

3)  Optimization  of  total  cost  subject  to  system  requirements  and  engineering  constraints 

The  given  optimization  problem  can  itself  have  many  parameters  which  can  be  in  one  of  the  below 
mentioned  categories. 

1)  Geometric  factors 

(i)  length, 

(ii)  cross-sectional  area, 

(iii)  number  and  location  of  electrodes, 

(iv)  electrode  connections 

(v)  Distribution  of  external  load  for  power  extraction 

(vi)  Extent  of  the  region  for  power  extraction 

(vii)  Connection  scheme  between  generator  and  accelerator 
(viii)  Efficiency  of  coupling  between  generator  and  accelerator 

2)  Physical  factors 

(i)  strength  of  B-field 

(ii)  electrical  conductivity  (depends  on  seeding  concentration,  or  beam  current,  beam  energy, 
beam  efficiency) 

Figure  below  shows  the  2-D  schematic  of  a  MHD  bypass  scheme  consisting  of  an  inlet,  generator, 
combustor,  and  accelerator.  The  non-dimensional  lengths,  the  reference  length,  typical  magnetic  fields 
and  electrical  conductivity  are  also  shown  in  the  figure. 


5.69  33  4.0 


Sigma  =  5  mho/m,  B  =  7  T 
Lref  =  0.6  m 
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The  following  figure  shows  an  optimization  problem  comprising  of  a  scheme  where  the  design 
parameters  are  the  coefficients  of  a  double  spline  for  both  the  generator  and  accelerator.  The  optimization 
problem  could  seek  to  maximize  the  extracted  power  in  a  generator  and  maximize  thrust  in  the 
accelerator.  The  some  possible  design  parameters  for  both  the  accelerator  and  generator  are  also  shown 
in  the  figure. 


Ig  =  length  of  generator 
l,  =  length  of  electrodes 
ag{x)  =  cross-sectional  area 

{aj.bpC^d^a^bj.CjjdjJg  —  co-efficients 
of  double  spline  for  generator 


la  =  length  of  accelerator 
i2  =  length  of  electrodes 
a/x)  =  cross-sectional  area 

{a1,b]>c,,d],a2,b2,C2,d2}a  =  co-efficients  of 
double  spline  for  accelerator 
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2.2  Shock-shock  interaction  and  stagnation  point  heating 

Avoidance  of  shock  wave  impingement  on  the  lower  cowl  lip  of  hypersonic  inlets  is  a  critical  problem  in 
the  design  of  hypersonic  flight  vehicles.  Shock/shock  interactions  can  result  in  an  amplification  of  the 
usual  stagnation  point  heat  transfer  rate  by  a  factor  of  30  or  more.  Edney1  reduced  shock/shock 
interactions  to  six  types  that  depend  on  the  orientation  of  the  shock  waves  to  the  body.  Only  types  III  and 
IV  involve  a  large  value  of  the  normalized  heat  transfer  rate,  q/q0,  where  q0  is  the  stagnation  point  rate 
without  any  interaction.  Type  III  and  IV  interactions  involve  jets  and  shear  layers  that  impinge  on  the 
surface  of  the  body  and  result  in  extreme  values  of  q/q0. 

One  can,  of  course,  design  the  inlet  ramps  so  that  shock  waves  do  not  impinge  upon  the  lower  cowl  lip  for 
the  entire  range  of  flight  mach  numbers,  however,  this  leads  to  excessive  spillage  drag.  Unfortunately,  use 
of  variable  geometry  to  control  shock  wave  position  is  often  not  feasible  in  hypersonic  flight. 

The  typical  optimization  would  involve  control  of  shock  position  via  MHD  to  minimize  spillage  drag 
(maximize  capture  ratio)  over  the  entire  flight  Mach  number  range,  without  allowing  the  shocks  to 
impinge  on  the  cowl  lip.  This  will  be  set  up  as  constraints  within  the  optimization  procedure.  The  shock 
impingement  will  be  forced  to  be  at  some  clearance  distance  away  from  the  cowl  lip.  A  simple 
implementation  of  this  strategy  would  be  to  increase  the  cowl  length  artificially  and  optimize  for  this 
configuration  so  as  to  maximize  the  mass  capture,  but  forcing  the  overall  configuration  to  cause 
intersection  of  the  shocks  at  the  tip  of  the  artificially  extended  cowl. 
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2.3  MHD-Based  Flow  control 


Turbulence,  transition,  drag  reduction. 


There  is  a  growing  resurgence  of  interest  in  the  use  of  MHD  to  influence  the  lift  and  drag  characteristics 
of  airfoils  and  hydrofoils.  The  control  of  turbulent  flows  and  transition  using  MHD  is  an  almost  classical 
MHD  problem,  which  has  reincarnated  multiply  in  recent  times.  Tsinober  has  provided  a  review  of  flow 
control  opportunities  available  here,  and  Henoch  et  al8  have  described  an  ongoing  experimental  program, 
which  holds  enormous  promise,  given  its  immediate  relevance  to  seawater  and  possible  extensions  to 
plasma  based  devices  in  aeronautical  engineering. 


Figure  30:  Lorentz  force  actuator  to  reduce  turbulent  skin  friction,  from  Ref.[8]  showing  side  and  top 
views.  Force  acts  in  the  spanwise  direction.  Electrodes  are  designed  to  reduce  spanwise  fringing  fields. 
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2.4  Optimal  MHD  Adjoint  Solution  via  the  Hamilton-Jacobi-Bellman  Equation 


A  generic  adjoint  based  optimization  procedure  is  outlined  for  an  Objective  Function,  Jc(U,<p,8 )  that 

depends  on  the  flow  variables  U,  electric  field  potential  (and  thereby  on  the  current)  <p,  and  the  design 
parameters,  6. 


Jc(U,<p,0)=  \p{U  ,<p,0)cKi  +  §M{U  ,<p,0)ds  (1) 

n  s 

For  example,  the  problem  could  be  the  maximization  of  the  mass  weighted  pressure  loss  for  a  hypersonic 

inlet  where  in  Jc(U,<p,0)  =  — — <^M{U ,(p,6)ds  and  M{u t<p,6)  =  puAy(P-Px)  or  even  the  total  mass  capture 

mP“  s 

at  the  inlet  to  avoid  any  spillage  drag  losses. 

6  is  the  ramp  angle  and  note,  in  this  case  along  the  stream-wise  direction  the  change  in  the  inlet  area  can 
be  represented  in  terms  of  the  ramp  angle  as,  Ay  —  (x  -  d)  tan  6,  where  d  is  the  stream-wise  location 
identifying  the  ramp  vertex.  The  vector  of  state  variables  (characterizing  the  flow)  are  given  by, 


'  p 

P 

P 

t/  = 

pu 

and  the  Flux  Vectors  in  Cartesian  coordinates,  F  = 

pu2+p 

G  = 

puv 

7 

pv 

puv 

pv  +p 

PE. 

(pE  +  p)u 

.( pE  +  p)v 

(2) 


~  +  R(u ,  6)  =  S(U,  cp,  6)  ^  ^  =  /({/,  e,  t )  where,  f(U,  0,  t)  =  -R{U,  0,  t)+  S(U,  <p,  9,  t ) 

ot _ dt_ _ 

(3) 

The  system  of  equations  governing  the  flow  (LHS)  without  the  source  term  is  purely  Hyperbolic,  while 
the  source  term  by  itself  is  computed  as  a  solution  to  an  Elliptic  PDE. 

We  adopt  an  approach  wherein,  the  flow  solution  is  first  computed,  setting  the  source  to  zero  and  then 
utilizing  the  flow  solution,  the  source  term  is  computed,  which  is  then  substituted  back  into  the  flow 
equations  for  the  solution  in  the  presence  of  a  non-zero  source  term  and  so  on. 


The  objective  function  however  depends  on  the  flow  variables,  source  term  states  and  the  design 
parameters. 


For  the  Flux  vectors  described  above,  the  Jacobians  with  respect  to  the  flow  variables  are  computed  to  be 

a F  ,  dG 

—  and  — . 
dU  dU 


Hamilton-Jacobi  Equation 

For  <?as  the  set  of  design  variables,  the  flow  fields  physics  is  represented  as  U(6)  where  U  is  assumed  to 
be  the  solution  of  (3).  Note  that  R(U,  6)  =  0,  is  the  steady  state  PDE. 

We  assume  that  the  computational  domain  and  the  primitive  variables  are  the  same  and  suppose  that  the 
cost  function  Jc  is  measured  b  (1). 
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This  optimal  cost  problem  is  viewed  in  the  sense  of  the  optimization  iteration  number  horizon,  where  a 
correspondence  between  the  iteration  number  and  the  time  horizon  in  conventional  optimal  control 
problems  is  hypothesized. 

The  justification  is  in  the  following  observation: 

Conventional  Optimal  Control  Problem: 

u(l)  =  ~K(x,t) 


is  the  feedback  law  that  is  sought  to  render  some  expected  cost  of  the  state  deviations  to  be  minimal.  Here 
u(t)  is  the  free  variable  (control)  that  is  obtained  as  a  function  of  the  states  of  the  dynamical  system. 

Flow  Control  Optimization  Problem:  The  desired  shape  change  is  expressed  as 
&k+ 1  =  &k  +  MU/c  ’<Pk)  ’ 


Where,  6k  is  vector  of  shape  parameters  at  the  k*  iteration  and  A (Uk ,  <pk )  is  a  function  of  the  states  of  the 
dynamical  system  (flow  solution). 


Thus  theoretically,  for  the  infinite  optimization  horizon  problem,  we  may  cast  this  problem  as  one  that 
seeks  the  optimal  control  strategy,  over  an  infinite  number  of  iterations. 


Interestingly,  these  two  problems  viewed  in  the  Discrete  Optimization  horizon  and  Discrete  Time  sense 
reduce  to  the  solution  of  the  Hamilton-Jacobi-Bellman  equation  that  is  summarized  below. 

For  the  system  described  in  (3),  with  the  steady  solution,  the  performance  index  can  be  summarized  as  in 

(1) 


The  value  function  is  then  defined  by, 


The  dynamic  programming  principle  states  that  for  every,  r  e  [1  N] 


.  o  r 

V'(U,<p,k)  =  —  Y  Jc  ( u ,  <p,  6)  +  V(r, U(r),  <p(r)) 
fry)  , 


From  the  above  we  obtain, 


VM  (U,  <p)  =  Vk  (£/,  cp) + //(¥,  Vy  V(U,  <p)) 
also 


dV * 
dt 


dv‘ 

5T 


/(%*)+ JC(¥,S) 


T  =  [[/r  <pf 


(4) 


V0(U,<p)  =  0 


The  Hamiltonian  H  is  given  as  //(T,  A)  =  [/c  (T,  3) + ir/(vF,  5)] 
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The  nonlinear  first  order  PDE  in  the  discrete  form  is  the  dynamic  programming  PDE  or  Hamilton-Jacobi- 
Bellman  equation.  Further,  the  Hamiltonian  is  concave  in  the  variable  A  (since  it  is  the  infimum  of  linear 
functions). 

Viscosity  Solutions 

The  terminology  for  the  viscosity  solution  comes  from  the  vanishing  viscosity  method,  which 
finds  solution  for  (4)  as  a  limit  of  V£  ->  F*  of  solutions  to 


~AV*  +F('¥,VC 


dv‘ 

5T 


■)  =  0 


(5) 


The  Laplacian  term  —  AVC  is  used  to  model  fluid  viscosity  (for  the  numerical  solution  of  the  HJB  alone). 


A  numerical  scheme  can  be  constructed  that  provides  the  optimal  shape  function  as  well  as  the 
Value  function  based  on  several  existent  solvers.  Posing  the  optimization  problem  as  an  optimal  control 
problem  this  way,  eliminates  the  need  to  solve  the  numerical  gradients.  Further,  in  case  of  the  Adjoint 
developments,  the  evaluation  of  the  Adjoints  is  not  trivial  especially  when  the  influence  of  shape 
functions  is  to  be  evaluated  on  a  far-field  flow  variable.  The  above  approach  alleviates  these  difficulties  at 
the  expense  of  the  iterative  computational  cost  (multiple  sub-problems).  The  relative  computational 
burden  should  lie  in  between  the  numerical  adjoint  solution  and  the  numerical  gradient  solutions. 

As  opposed  to  the  traditional  and  adjoint  methods,  where  the  focus  is  on  trying  to  compute  the 
gradient  of  the  cost  function  with  respect  to  the  shape  parameters,  the  focus  of  the  solution  to  the  HJB 
equation  is  to  construct  an  approximation  to  a  sequence  of  shape-changes  that  will  lead  to  the  minimum 
cost. 
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This  research  considers  the  effects  of  magnetohydrodynamic  source  terms  on  the 
optimization  of  a  scramjet  inlet.  A  representative  case  is  found  by  determining  the 
optimal  inlet  configuration  for  the  inviscid  case  assuming  no  MHD  source.  This  is  a 
simplified  look  into  the  optimization  of  the  two  ramp  inlet  say  to  a  supersonic  scramjet 
engine.  The  parameter  being  optimized  is  the  mass  weighted  pressure  loss  at  the  throat 
of  the  inlet  such  that  the  pressure  drop  through  the  inlet  is  minimized.  The  optimal  inlet 
configuration  is  then  determined  allowing  for  the  presence  of  a  charged  flow  through  a 
magnetic  field.  This  comparison  will  demonstrate  the  effects  of  the  MHD  source  on  the 
flow. 
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CHAPTER  3 


OPTIMIZATION 

Optimization  is  the  process  that  finds  the  best  or  optimal  solution  for  a  given 
situation.  There  are  three  basic  ingredients  that  make  up  an  optimization  problem: 

1.  An  objective  function  f(x )  -  That  which  we  want  to  minimize  or  maximize. 

2.  A  set  of  variables  or  unknowns  x  =  {aq,  x2, . . .  xn}  -  The  parameters  that  will  affect 
the  value  of  the  objective  function. 

3.  A  set  of  constraints  -  Limits  set  on  the  variables  in  any  of  the  following  forms: 

•  equality  — >  Gi(x)  =  0,  (i  =  1, . . . ,  me ) 

•  inequality  — >  Gi(x)  <0,  (i  =  me  +  1, . . . ,  m) 

•  parametric  bounds  — >  ximueT,  xupper 

The  optimization  problem  can  then  be  stated  mathematically  as: 

min  f{x) 
xeJS"  ' 

subject  to 

Gi(x)  =  0  (i  =  l,...,me) 

G,(x )  <  0  (i  =  me  +  1, . . .  ,m) 

Slower  —  —  %upper 

The  type  of  optimization  routine  chosen  to  solve  a  particular  problem  is  dependant 
on  the  size  and  characteristics  of  that  problem.  The  size  is  influenced  by  the  number  of 
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design  variables  and  constraints,  and  the  characteristics  describe  the  linear  or  nonlinear 
nature  of  the  problem,  this  applies  to  the  objective  function  as  well  as  the  constraints. 

In  our  case  the  design  variables  are  the  two  ramp  angles,  x  =  {01,0:2}.  The 
objective  function,  F0bj(x)  we  are  trying  to  optimize  is  the  mass  weighted  total  pressure 
loss: 


Fob] 


Y2  Pi,jui,j^yi,j(F>i,j  Poo) 

rhPoc 


(3.1) 


where  the  i,  j  indices  refer  to  the  mesh  points  on  the  exit  plane  of  the  inlet  and  Phj 
and  Poo  are  the  total  pressures  at  the  exit  mesh  points  and  at  the  entrance  to  the  inlet 
respectively.  And  the  constraints  are  the  flow  equations,  solved  by  use  of  the  flow  solver 
as  well  as  the  condition  that  a\  <  a2-  We  also  set  some  upper  and  lower  limits  on  the 
angles  to  ensure  a  well  behaved  result. 

Our  problem  is  therefore  a  constrained  nonlinear  multivariable  problem  in  which 
we  are  going  to  implement  a  gradient  based  method.  This  type  of  method  is  generally 
the  most  efficient  when  the  function  to  be  minimized  is  continuous  in  it’s  first  derivative. 
For  functions  that  are  highly  nonlinear  or  have  discontinuities,  a  search  method  that 
utilizes  only  function  evaluations  is  more  suitable.  Where  as  higher  order  methods  such 
as  Newton’s  method,  may  be  feasible  for  functions  where  the  second  order  information 
is  easily  calculated. 


3.1  BFGS  Optimization 

As  mentioned,  Newton’s  method  is  a  higher  order  method  that  requires  the  calcu¬ 
lation  of  Hessian,  H  at  each  step.  This  is  a  very  computationally  expensive  procedure 
and  can  at  times  be  very  difficult.  We  there  for  look  to  a  set  of  methods  known  as 
quasi-Newton  or  variable  metric  methods.  Here  an  approximate  Hessian  matrix  is  built 
up  gradually  instead  of  calculating  it  at  every  single  point.  The  build  up  of  the  Hessian  is 
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done  using  gradient  information  from  some  or  all  of  the  previous  iterations  to  formulate 
a  quadratic  model  problem  of  the  form: 

min  \xtHx  +  cTx  +  b 

x  2 

where  the  Hessian  H  is  a  positive  definite  symmetric  matrix,  c  is  a  constant  vector  and  b 
is  a  constant.  The  optimal  solution  for  this  problem  occurs  when  the  partial  derivatives 
of  /  go  to  zero, 

V/(x*)  =  Hx*  +  c  =  0 

where  the  optimal  solution  x*  can  be  written  as: 

x*  =  H~'c 

Quasi-Newton  methods  use  the  observed  behavior  of  f(x )  and  V  f(x)  to  build  up 
curvature  information  to  make  an  approximation  to  H.  The  Broyden-Fletcher-Goldfarb- 
Shanno  (BFGS)  method  is  just  one  technique  for  the  calculating  the  Hessian  update.  It 
is  considered  to  be  one  of  the  most  effective  general  purpose  methods  [26]. 

The  BFGS  optimization  formula  involves  calculating  a  search  direction  in  order  to 

determine  the  descent  direction  that  will  locate  the  minimum. 

\ 

dk  =  -H;1  ■  V/(x*)  (3.2) 

where  f(x)  if  the  objective  function  and  k  denotes  the  iteration  number.  We  start  with 
an  initial  design  vector  x0  and  an  approximate  Hessian,  H0  which  can  be  set  to  any 
symmetric  positive  definite  matrix,  i.e.  identity  matrix. 
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A  one-dimentional  search  is  performed  in  the  search  direction  to  determine  the 


distance,  a*  to  be  traversed  in  that  direction.  The  design  is  then  updated  as: 


Xk+ 1  —  Xjc  T  Oi  dk 


(3.3) 


The  Hessian  matrix  is  then  updated  as: 


Hk+ i  —  Hk  + 


gfcgfc 

€sk 


HTkslskHk 

sjHkSk 


where 


(3.4) 


Sjc  —  %k+ 1  %k 

qk  =  V/(x*+i)  -  V/(xfc) 

Gradient  information  is  either  supplied  through  analytically  calculated  gradients, 
or  derived  by  partial  derivatives  using  a  numerical  differentiation  method  via  finite  dif¬ 
ferences.  This  involves  perturbing  each  of  the  design  variables  in  turn  and  calculating 
the  rate  of  change  in  the  objective  function. 

The  BFGS  method  implemented  here  is  in  the  form  of  the  ‘fmincon’  function  in 
MATLAB®  \  ‘fmincon’  is  designed  to  find  a  minimum  of  a  constrained  nonlinear  multi- 
variable  function  [26].  It  uses  a  Sequential  Quadratic  Programming  (SQP)  method  where 
an  estimate  of  the  Hessian  of  the  Lagrangian  is  updated  at  each  iteration  using  the  BFGS 
formula.  The  Lagrangian  is  the  cost  function  augmented  with  the  constraints: 

L(x,  A)  =  F(x)  -  \TC(x) 


'MATLAB/SIMULINK  is  the  trademark  of  The  Math  Works  Inc;  Natick,  MA,USA 
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The  function  ‘fmincon’  was  found  to  be  very  sensitive  to  the  step  size,  as  well  as 
the  scaling  of  the  objective  function  itself.  Several  combinations  of  the  above  mentioned 
variables  were  tried  in  order  to  obtain  the  best  results  and  it  was  found  that  a  scaling  of 
ln(  1  +  Fobj 2)  provided  good  convergence  properties. 

In  order  to  ensure  that  we  are  still  minimizing  the  correct  function  we  state  the 
following  lemma: 

LEMMA  1  (Modified  Objective) 

min/n(l  +  M L)  =  min  Fobj 

X  J  X 


Proof: 

To  find  an  extremum  of  a  function  we  take  the  first  derivative  and  set  it  equal  to  zero: 

of 

This  equation  is  only  satisfied  when  Fobj  —  0,  which  is  a  trivial  solution,  or  when  =  0, 
which  is  the  same  extremum  as  the  original  function,  minxF0y. 

Now  we  wish  to  prove  that  this  extremum  is  in  fact  the  same  extremum  as  that 
of  the  original  function,  we  therefore  take  the  second  derivative  and  confirm  that  it  is  of 
the  same  sign: 


/  2 F0bj  \  dFpbj  f  2 Fobj  \  d2Fobj 

dx  \l  +  Fobj )  dx  \l  +  F2obj J  dx2 

(  2Fobj  \  d2Fobj 

W  +  f%) 


33 

since  Fabj  is  a  positive  value,  the  result  is  simply  a  scaling  of  the  second  derivative  of 
the  original  function,  the  sign  is  unchanged.  This  proves  that  extremum  of  modified 


objective  function  is  of  the  same  type  as  that  of  F^j,  therefore: 

in  in  ln(  1  +  F„bj)  =  min  Fobj 


3.2  Implementation 

The  implementation  of  this  optimization  routine  is  laid  out  in  figure  3.1.  As  is 
shown  in  the  flow  chart  the  flow  solution  is  calculated  n  +  1  times  for  each  update  of 
the  Hessian,  where  n  is  the  number  of  design  parameters.  The  values  of  the  objective 
function  are  stored  from  iteration  to  iteration  for  determination  of  the  minimum  and 
calculation  of  the  gradients. 

3.3  Results 

The  final  optimized  results  of  =  3.904  deg  and  a2  =  7.265  deg  are  consistent 
with  Korte’s  results  of  a\  =  4.263  and  a2  =  7.621.  The  initial  and  final  configurations 
are  summarized  in  table  3.1.  The  initial  condition  was  the  shock  cancelled  case  discussed 
earlier. 

Fig.  (3.2)  shows  the  design  history  of  the  optimization  routine  at  each  iteration  of 
the  inside  loop  in  Fig.  (3.1).  The  outside  loop  iterations  or  number  of  times  that  the 
Hessian  was  updated,  is  equal  to  the  total  number  of  flow  solver  iterations  divided  by 
the  number  of  design  variables  plus  1,  (^j).  From  Fig.  (3.1)  we  can  tell  that  there 
were  approximately  93  iteration  of  the  flow  solver  and  2  design  parameters,  giving  us 
approximately  31  update  of  the  Hessian.  The  Fabj  values  plotted  are  scaled  as  per  the 
scaling  function  mentioned  earlier  so  as  to  improve  the  convergence  properties.  Fig.  (3.3) 
is  a  plot  of  the  Mach  number  contours  for  the  final  optimized  configuration. 


Figure  3.1  Flowchart  of  optimization  routine. 


Table  3.1  Summary  ol 


Initial  condition 
Optimized  value 


initial  and  final  inlet  conditions. 
QU  (deg)  a2  (deg)  Fobj 
2.882  9.342  0.2445 

3.904  7.265  0.1718 


Angles  ci^,  a2,  and  Objective  Function  Vs  Optimizer  Iteration 


t - 1 - 1 - 1 - 1 - 1 - -“i - r 


_J _ 1 _ 1 _ I _ I _ I _ I _ L_ 

20  30  40  50  60  70  80  90 

iteration 


Figure  3.2  Optimization  design  history  of  ramp  angles. 


CHAPTER  4 


INVISCID  EULER  EQUATIONS  WITH  MHD  SOURCE 
4.1  MHD  source  term  modelling 

The  effects  of  a  charged  flow  through  the  inlet  with  an  applied  electromagnetic 
field  can  be  modelled  with  the  addition  of  appropriate  electrodynamic  terms  [7].  Let  us 
first  consider  a  charged  flow  with  conductivity,  a  and  velocity  V;  if  we  apply  a  magnetic 
field,  B  to  this  flow  the  resultant  current  density  is  given  as  Jb  =  crV  x  B.  Whereas  if 
an  electric  field,  E  is  applied  to  the  same  flow  the  current  density  is  given  as  Je  =  crE. 
The  contributions  to  the  electric  field  may  be  due  to  an  applied  electrostatic  field  Es, 
due  to  the  mutual  repulsion  or  attraction  of  electric  charges,  or  an  induced  electric  field, 
Ej,  due  to  the  flow  of  charged  particles  through  a  magnetic  field,  such  that  E  =  Es  +  E;. 
Thus  the  total  current  density  in  a  flow  with  an  applied  electromagnetic  field  (neglecting 
Hall  effect)  is  given  by  the  generalized  Ohm’s  law: 

J  =  cr(Es  +  Ej  +  V  x  B)  (4.1) 

For  the  electrostatic  field  we  know  from  Coulomb’s  law  that  Es  is  irrotational  and  from 
Gauss’s  law  that  the  divergence  is  fixed: 

V  x  Es  =  0 
V  ■  Es  =  ^ 

£o 
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where  pe  is  the  charge  density  and  e0  is  the  permittivity  of  free  space.  The  induced  electric 
field  on  the  other  hand  has  zero  divergence  and  a  finite  curl  governed  by  Faraday’s  law: 


V  •  E; 


V  x  E; 


0, 

_dB 

dt 


We  therefore  have  a  total  electric  field  E  =  Es  +  Ej,  uniquely  determined  by: 

V  •  E  =  — ,  (Gauss’s  law)  (4.2) 

£o 

SB 

V  x  E  =  —  — ,  (Faraday’s  law)  (4.3) 

and  J  =  er(E  +  V  x  B)  (Generalized  Ohm’s  law)  (4.4) 


MHD  flows  are  characterized  as  having  a  very  low  (<  1)  magnetic  Reynolds  number,^ 


Rm  —  fJ'oO'gV'L' 

where  is  the  permeability  of  free  space,  ag  is  the  gas  conductivity,  u  is  the  gas  velocity 
and  L  is  a  characteristic  length. (High  magnetic  Reynolds  numbers  (>  l)are  characteristic 
of  fusion  research  and  astrophysical  phenomena).  This  tells  us  that  the  conductivity  in 
MHD  flows  is  very  low,  therefore  the  currents  and  hence  the  induced  electric  field  is  also 
very  small.  This  allows  us  to  assume  B  to  be  constant,  therefore 

V  xE=~«0 

at 
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As  a  consequence  of  the  above  we  may  introduce  a  scalar  electrostatic  potential  p  such 
that  E  =  ~V<p  and  VV  —  const.  We  can  then  replace  the  electric  field  term  in  the 
generalized  Ohm’s  law,  Eq.  (4.4)  with  the  electrostatic  potential  to  get 


J  =0{~V<p  +  V  x  B) 


(4.5) 


Since  V  •  (V  x  E)  =  0,  Faraday’s  law,  Eq.  (4.3)  implies  that  the  magnetic  field,  B  is 
solenoidal;  V  ■  =  0,  therefore  V  •  B  =  0. 

The  conservation  of  charge  is  given  as, 


V  J  =  - 


dpe 

dt 


this  equation  simply  states  that  the  rate  at  which  the  charge  is  decreasing  in  a  small 
volume  must  equal  the  rate  at  which  the  charge  flows  out  across  the  surface  of  that 
volume.  For  a  stationary  conductor  pe  is  zero,  and  when  there  is  motion  it  turns  out  that 
pe  is  very  small  and  too  low  to  produce  any  significant  force,  pe E  and  therefore  may  be 
neglected.  Based  on  these  observations  we  conclude  that 


V- J  =  0 


therefore,  J  is  also  solenoidal.  Taking  the  divergence  of  Eq.  (4.5).  we  can  further  simplify 
the  generalized  Ohm’s  law  to  get  an  expression  for  the  electric  potential 


V- J  =  0 


=»  V  •  <r(Vy?)  =  V  •  (<jV  x  B) 
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If  we  integrate  both  sides  of  this  equation  over  a  volume  dV  and  use  Gauss’s  theorem  to 
convert  the  volume  to  a  surface  integral  the  equation  becomes, 


f  V  •  cr(Vip)dV  =  f  V  •  (ffV  x  B )dV 
Jv  Jv 


J  aV(pds  =  J 


crV  x  Bds 


We  can  now  discretise  this  formula  for  use  in  our  flow  solver  to  calculate  the  MHD 
portion  of  the  flow  simply  by  summing  up  the  contributions  from  each  grid  cell  surface. 
For  a  given  magnetic  field  and  flow  conductivity  we  can  calculate  the  electric  potential 
as  follows: 


E  "1  ^  =  E  */<V  -  B>A* 

faces  '  '  faces 


E,**.  <7(V  x  B)As  -  As 

Vp  = - V - SAi - 

faces  d 


This  allows  us  to  then  calculate  the  current  density  directly  from  Eq.  (4.5),  and  it’s 
impact  on  the  flow  equations  is  demonstrated  in  the  following  section. 


4.2  Euler  equations  with  MHD  effects 

As  derived  in  chapter  2  the  two-dimensional,  nonlinear,  Euler  equations  in  conser¬ 
vative  form  for  a  compressible,  inviscid  flow  are  given  as 


dU  dF  dG 
dt  +  dx  dy 


(4.6) 


With  the  application  of  a  magnetic  field  to  a  charged  flow  the  body  forces  and  volumetric 
heating  effects  are  no  longer  negligible.  The  body  force  term  known  as  the  Lorentz  force 
are  given  by  the  vector  J  x  B,  and  volumetric  heating  known  as  Joule  heating  is  given 
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*2 

as,  As  seen  earlier  in  the  development  of  the  Euler  equations  the  contribution  to  the 

momentum  equation  is  strictly  due  to  the  Lorentz  forces  where  the  contribution  to  the 

-2 

energy  equation  is  the  total  rate  of  energy  addition,  J  •  E  —  +  V  •  (j  x  B)  due  to  both 

volumetric  heating  and  work  done  by  the  Lorentz  forces. 

In  Eq.  (4.6)  U,  F  and  G  are  as  defined  in  the  previous  section  and  the  source  term, 
S  is  given  as 

0 

(J  x  B)x 
(J  x  B)y 
(J  '  E) 

The  left  hand  side  of  Eq.  (4.6)  is  unchanged  from  before,  and  as  stated  earlier  this 
equation  is  hyperbolic  for  Mach  numbers  greater  than  one.  The  right  hand  side  however 
is  unconditionally  elliptic  for  smooth  variations  of  material  properties,  we  therefore  need 
to  implement  a  Poisson  solver  to  obtain  the  source  terms. 

4.2.1  Magnetic  field  perpendicular  to  the  x-y  plane 

For  a  magnetic  field  applied  to  a  conducting  fluid,  there  is  an  induced  electric  field. 
Knowing  this  fact  we  must  allow  for  the  influence  of  the  Hall  effect  in  our  equations,  such 
that  Ohm’s  law  must  be  rewritten  as  follows 

J  =  <t(-Vv?  IVxB)-^(JxB) 

B 


(4.7) 
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the  term  uit  is  called  the  Hall  parameter  where  r  is  the  mean  free  time  between  collisions, 
and  u  is  the  cyclotron  frequency  which  for  an  electron  with  charge  e  and  mass  me,  moving 
through  a  magnetic  field  of  strength  B,  is  equal  to 


(4.8) 


With  the  magnetic  field  implemented  along  the  z-axis  and  utilizing  the  property  of 
J  being  divergence  free,  we  can  convert  Eq.  (4.7)  into  a  differential  equation  in  terms  of 
the  electric  potential.  Let  us  first  set  the  following  notation: 


V  =  (u,v,0),  B  =  (0, 0,  B)  (3  =  lot,  a  =  i  + 


Then  we  can  break  down  Eq.  .(4.7)  into  it’s  components: 


J X 

a{~ <px  +  Bv)  -  pjy 

Jy 

= 

cr(~ipy  -  Bu )  +  pjx 

JZ 

0 

Since  we  do  not  apply  an  external  electrostatic  field  the  only  electric  field  present  is  due 
to  the  applied  magnetic  field  and  is  restricted  to  the  x-y  plane  thus  the  z-component  of 
the  electric  potential,  tpz  is  zero.  The  current  density  components  are  therefore  given  by 


Jx 

<px  —  (3tpy  —  Bv  —  /3Bu 

Jy 

=  a 

^Py  +  (3ipx  +  Bu-  f3Bv 

and  if  we  now  enforce  the  divergence  free  condition  on  the  current  density,  V  •  J  =  0 


d_ 

dx 


grouping  like  terms 
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d_ 

dx 


9<P  nfy 

a  q - Oij3 — 

ax  dy 


d_ 

dy 


r\  r\ 

=  ■£- a(vB  +  (3uB)  -  a(uB  -  pvB) 
ox  oy 


if  we  assume  constant  properties,  a  and  (5  can  be  pulled  out  of  the  partial  derivatives 


dV  _  V 

8x 2  ^dxdy 


+  aft 


dtp2 

dxdy 


vV 


r\  r\ 

—a(vB  +  puB)  —  —a(uB  —  fivB) 
dx  dy 

r\  r\ 

^  (vB  +  puB)  -  —  (uB  -  (4.9) 


Eq.  (4.9)  is  the  final  form  of  the  Poisson  equation  to  be  solved  for  this  case. 


4.2.2  Magnetic  field  in  the  x-y  plane 

If  however,  we  implement  a  magnetic  field  in  the  x  —  y  plane  and  consider  an  ideally 
sectioned  Faraday  MHD  generator  such  that  the  Hall  effect  is  neutralized,  the  resultant 
current  density  is  only  in  the  z  direction: 


V  =  (u,u,0),  B  =  {Bx,By,0),  J  =  (0, 0,  J) 
J  =  cr(- Vv?  +  V  x  B) 


1 _ 

0 

Jy 

=  <7 

0 

1 - 

uBy  -  vBx  -  ipz 

For  this  research  we  are  setting  ipz  to  zero  which  corresponds  to  a  short  circuit  of  the 
electric  field. 
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The  source  terms  are  then  trivially  found  as  follows  without  the  need  for  a  poisson 


solver: 


0 

CTBy(i/Bx  —  uBy) 
aBx(uBy  —  t>Bx) 

J  E 


4.3  Results 

4.3.1  Magnetic  field  perpendicular  to  the  x-y  plane 

To  simulate  an  ionized  flow  we  set  the  conductivity  in  a  small  pillbox  area  in  the 
center  of  the  inlet  where  the  flow  would  most  likely  become  ionized,  to  some  chosen 
value,  and  the  conductivity  outside  the  box  is  set  to  le~6.  For  this  case  we  oriented  the 
magnetic  field  along  the  z-axis  or  out  of  the  page  with  a  uniform  magnetic  field  strength 
of  1.0  tesla.  As  you  can  see  in  Fig.  (4.1)  the  currents  are  indeed  generated  but  only 
in  the  zone  of  simulated  ionization,  this  is  expected  since  the  conductivity  is  negligible 
outside  that  zone.  We  used  a  full  inlet  in  this  case  due  to  the  non-symmetric  nature  of 
the  problem.  Fig.  (4.2)  demonstrates  the  effects  that  this  current  flow  has  on  the  Mach 
contours,  the  flow  decelerates  earlier  and  the  pressure  recovery  at  the  inlet  is  decreased 
from  the  non-MHD  case. 

4.3.2  Magnetic  field  in  the  x-y  plane 

For  a  magnetic  field  implemented  in  the  x-y  plane  there  are  no  currents  to  be  seen 
since  they  are  generated  perpendicular  to  the  x-y  plane.  We  ran  a  number  of  simulations 
with  varying  conductivity  and  magnetic  field  angle  and  these  results  are  shown  in  Fig. 
(4.3)  for  two  different  magnetic  field  strengths. 
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As  can  be  seen  in  Fig.  (4.3)  the  application  of  a  magnetic  field  to  the  charged 
flow  always  results  in  a  decrease  in  the  pressure  recovery  (or  an  increase  in  the  pressure 
loss).  This  is  consistent  with  the  theory  in  that  we  are  not  adding  any  energy  to  the 
flow  and  the  resultant  joule  heating  is  only  a  detriment  to  the  pressure  recovery.  We 
therefore  adjust  the  problem  to  consider  a  cowl  style  inlet  where  we  try  to  optimize  the 
mass  capture.  The  optimal  situation  for  the  cowl  configuration  can  be  seen  in  Fig.  (4.4). 


Cowl 


Figure  4.4  Optimal  cowl  configuration  with  shocks  converging  on  the  cowl  lip. 

To  determine  our  optimal  configuration  we  ran  the  optimizer  with  a  mass  capture 
objective  function  given  as: 

F0bj  —  ^  )  pjjUjj&yij  (4.10) 

The  results  of  which  are  shown  in  Fig.  (4.5)  where  a\  =  2.2  deg,  <22  =  8.9  deg  and  the 
mass  capture  equal  to  6.1965  kg/s. 

During  off  nominal  flight  conditions  where  the  flow  Mach  number  is  less  than  the 
design  value,  the  shocks  will  move  ahead  of  the  cowl  lip  and  some  of  the  compressed  air 
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Figure  4.5  Mach  contours  of  optimized  cowl  inlet. 


will  escape  the  inlet  resulting  in  ‘spillage’  and  a  decrease  in  the  mass  capture.  In  flight 
conditions  where  the  flow  Mach  number  is  greater  then  the  design  value,  the  shocks  move 
into  the  inlet  causing  multiple  reflected  shocks,  loss  of  total  pressure,  possible  boundary 
layer  separation  and  engine  unstart  [4].  Fig.  (4.6)  demonstrates  these  conditions. 

To  simulate  a  less  than  design  Mach  number  flow  we  adjust  the  ramp  angles  to 
ai  =  3-0  deg  and  a.2  —  9-0  deg  which  results  in  a  mass  capture  of  5.8757  kg/s.  See  Fig. 
(4.7)  for  the  Mach  contours  in  this  case. 
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Inlet  flow  Mach  Number  Contours 


Figure  4.7  Mach  contours  for  off  nominal  cowl  inlet. 


We  now  investigate  the  ability  of  an  applied  magnetic  field  to  direct  the  flow  back 
to  the  optimal  mass  capture  configuration.  Two  different  scenarios  are  considered;  one  is 
a  moving  e-beam  type  ionization  method  [3]  such  that  the  magnetic  field  and  ionization 
region  are  movable  and  coincident,  the  second  is  a  moving  magnetic  field  but  stationary 
ionized  region  [4],  Fig.  (4.8)  demonstrates  each  scenario.  In  the  case  of  the  e-beam 
ionization  the  ionized  region  is  enclosed  by  within  lines  that  are  parallel  to  the  magnetic 
field  lines:  yo  +  (x  —  xo  —  Ax)  tan(a)  <  y(x)  <yo  +  (x  —  Xo  +  Ax)  tan(a),  the  line  y(x)  < 
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radius  of  inlet,  and  the  body.  The  parameters  x0 ,  y0  and  Ax  determine  the  position  and 
width  of  the  ionized  region. 

Figs.  (4.9)  through  (4.12)  show  the  results  of  this  study.  It  is  clear  from  these 
results  that  the  larger  the  magnetic  field  strength  and  the  larger  the  conductivity  the 
greater  the  influence  on  the  flow.  It  is  interesting  to  note  that  the  stationary  ionized  zone 
has  a  much  greater  ability  to  manipulate  the  flow  than  the  coincident  ionized  region.  For 
the  largest  values  of  B  and  a  (line  corresponding  to  B  =  0.5  in  Fig.  (4.12)),  we  can  see 
that  the  mass  capture  is  indeed  approaching  that  of  the  optimal  situation. 

With  the  results  from  this  broad  parametric  study  we  decided  to  set  up  the  prob¬ 
lem  in  the  optimization  routine  and  compare  the  results.  This  time  we  set  the  objective 
function  as  the  mass  capture  Eq.  (4.10),  and  the  design  variable  was  the  magnetic 
field/electron  beam  angle  alone.  Aware  of  the  sensitivities  in  the  optimization  routine 
we  decided  to  try  several  different  initial  conditions  to  see  how  well  the  results  con¬ 
verged.  We  decided  on  initial  conditions  of  the  implementation  angle,  Qsi  above  and 
below  the  approximate  optimal  values  given  by  the  parametric  study.  Tables  (4.1)  and 
(4.2)  summarize  the  results. 

Table  4.1  Table  of  optimized  magnetic  field  angle  for  coincident  ionized  zone. 


0Bi  below  optimal 

0Bi  above  optimal 

Initial  condition  (deg) 

100 

170 

Optimized  value  (deg) 

141 

165 

Table  4.2  Table  of  optimized  magnetic  field  angle  for  stationary  ionized  zone. 


&Bi  below  optimal 

0Bi  above  optimal 

Initial  condition  (deg) 

80 

Optimized  value  (deg) 

108 

111 

Magnetic  field  angle  (deg) 
(b)  Stationary  ionization 


Figure  4.9  Mass  capture  for  inlet  with  conductivity  a  =  0.5  mho/m.  Results  are  shown 
for  (a)  coincident  ionization  and  (b)  stationary  ionization. 
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58 

Comparing  the  results  of  the  optimizer  with  those  of  the  parametric  study  for  the 
stationary  ionized  zone  we  can  see  a  very  nice  correlation.  The  spread  in  the  optimized 
results  is  not  large  and  is  consistent  with  the  parametric  study.  However  in  comparing 
the  optimizer  results  with  those  of  the  parametric  study  in  the  case  of  the  coincident 
ionization  zone  we  see  a  bit  of  discrepancy.  There  is  a  large  spread  in  the  values  given  by 
the  optimizer  depending  on  whether  we  began  the  optimization  above  or  below  the  peak 
value  shown  in  the  parametric  study.  In  order  to  better  understand  this  result  we  did 
another  parametric  study  around  the  peak  value.  We  limited  the  study  to  the  B  —  0.5 
tesla  and  a  =  2  mho/m  case  and  varied  the  magnetic  field  angle  from  120  to  180  degrees. 
As  can  be  seen  in  Fig.  (4.13)  the  results  of  this  study  show  that  the  curve  is  not  very 
smooth  and  this  might  explain  the  large  range  in  values  given  by  the  optimizer  for  this 
case. 

To  round  out  the  study  we  decided  to  try  to  optimize  both  the  geometry  and  the 
magnet  field  angle.  This  case  may  seem  redundant  in  that  it  would  require  both  actuating 
surfaces  as  well  as  an  ionization  method  and  a  way  to  generate  the  magnetic  field,  but  it 
is  also  represents  the  possibility  of  fine  tuning  the  magnetic  field.  Again  the  off  nominal 
case  of  cti  =  3.0  deg  and  a2  =  9.0  deg  was  used  with  the  initial  magnetic  field  angle  of 
70  deg.  We  limited  this  study  to  the  stationary  ionization  zone  with  a  conductivity  of  2 
mho/m  and  a  magnetic  field  strength  of  0.5  tesla,  the  results  of  which  are  shown  in  Fig. 
(4.14). 

The  final  configuration  for  the  geometry  is  very  close  to  that  of  the  optimal  mass 
capture  with  no  MHD  source  term  present,  and  the  final  magnetic  field  angle  is  consistent 
with  that  of  the  previous  case.  The  results  are  summarized  in  table  4.3. 


Magnetic  field  angle  (deg) 


Figure  4.13  Zoom  in  of  peak  values  for  the  moving  ionization  zone. 


Table  4.3  Table  of  optimized  geometry  and  magnetic  field  angles  for  stationary  ionized 

zone. 


c*i  (deg) 

a2  (deg) 

8b  (deg) 

Initial  condition 

3.0 

9.0 

70.0 

Optimized  value 

2.14 

9.11 

110.5 
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